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The phase behavior of the Baxter adhesive hard sphere fluid has been determined using specialized 
Monte Carlo simulations. We give a detailed account of the techniques used and present data for the 
fluid-fluid coexistence curve as well as parametrized fits for the supercritical equation of state and 
the percolation threshold. These properties are compared with the existing results of Percus-Yevick 
theory for this system. 
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I. INTRODUCTION 

For a pure substance to be capable of exhibiting equi- 
librium between two fluid phases — a gas and a liquid — 
the interaction between its particles must include an at- 
tractive component. Possibly the simplest model to in- 
corporate attraction is the square-well potential, where 
particles present each other with a hard core and a uni- 
form region of favorable potential energy up to an abrupt 
distance threshold beyond which there is no interaction. 
The square- well model has two intrinsic length scales: the 
hard-core diameter and the width of the attractive well. 
The ratio of these lengths determines the phase behavior. 

In 1968, Baxter introduced a model fluid containing 
attraction but having only one length scale, the hard- 
core diameter.^ In this model, the attraction takes the 
form of an adhesive interaction when particles are actu- 
ally in contact, and therefore has an effectively infinites- 
imal range. Baxter produced an analytic equation of 
state for this "adhesive hard sphere" model using the 
Percus-Yevick framework and compressibility equation, 
showing that the model possesses a critical point below 
which gas-liquid equilibrium may be observed. A sec- 
ond Percus-Yevick equation of state — derived from the 
energy equation — soon followed)^ with a strikingly differ- 
ent critical point from the compressibility-equation pre- 
diction. 

The Percus-Yevick framework also furnishes the clus- 
ter size distribution^ The surface adhesion provides an 
unambiguous definition of clusters as connected sets of 
mutually touching particles. At sufficiently low density, 
the fluid consists only of monomers and finite clusters, 
but above a temperature-dependent density threshold, 
clusters tend to aggregate into amorphous macroscopic 
structures. The system then contains sample-spanning 
infinite clusters and is said to percolate. The percola- 
tion threshold, at which the divergence of the cluster size 
occurs, has been located within the Percus-Yevick ap- 
proximation by Chiew and Glandtn^ 

Baxter regarded the adhesive hard sphere model as a 
highly idealized manifestation of a general potential with 
a repulsive core and an attractive tail. However, it has 



since been suggested^ that it is quite a reasonable descrip- 
tion of colloidal systems, where the particles interact on 
length scales much smaller than their own size. The gas 
and liquid are then more appropriately regarded as low- 
and high-density colloidal fluid phases. Percolation can 
be an observable phenomenon in colloidal systems, and 
may sometimes be detected directly, such as by a sud- 
den increase in electrical conductivity across the sample>^ 
The Baxter model is particularly appealing in the light of 
recent theoretical^ and experimental^ work that predicts 
and confirms the existence of a re-entrant glass-liquid- 
glass transition in colloidal systems with short-range at- 
tractive forces. For very short-range attraction there is 
also a glass-glass transition between a caged "repulsive" 
glass and an energetic "attractive" glass. 



A major obstacle in the application of the adhe- 
sive hard sphere model to the description of experi- 
mental data is that the model's phase behavior is only 
known through the approximate and conflicting the- 
oretical results of Percus-Yevick theory. Computer 
simulationsSiSiiSiii have provided some information on 
percolation, but have not tackled the issue of phase co- 
existence. The purpose of the present contribution is to 
provide a better knowledge of the adhesive hard sphere 
phase diagram using computer simulation. We describe 
in detail the special Monte Carlo techniques employed 
(Section mil) and supply numerical results for the fluid- 
fluid coexistence curve, the percolation threshold, and 
the supercritical equation of state (Section Hvl) . We be- 
gin, in Sectional with a definition of the model itself. 
It will be seen that this definition involves taking a limit 
that introduces a singularity into the potential energy 
function — an alarming prospect at first sight. The model 
in its strictest interpretation is indeed pathological,^^ but 
the problem is rather more subtle than it first appears, 
and in fact starts to set in before the full limit is taken. 
A proper discussion of this point is deferred to Section 
nil Dl where the difficulties can be put into the context 
of the simulation methodology. 
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II. THE ADHESIVE HARD SPHERE MODEL 

The adhesive hard sphere model is defined by the pair 
potential functioni 
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in the limit d — > cr, where kT is the thermal energy and 
r the particle separation. Eq. ^ describes a square well 
with hard-core diameter cr, width d — cr, and depth con- 
trolled by the parameter r. 

Equation couples the depth of the well to its width 
in such a way that the total Boltzmann weight of bound 
configurations remains finite as the well width is reduced 
to zero: 

lim Te-^W/'^^Wdr = lim -^^{d^ + ad + a^) 
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Hence, in the limit d ^ cr, which corresponds to an 
infinitesimally narrow but infinitely deep well, one can 
write 
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The Heaviside step function, 8, accounts for the hard- 
core repulsion, while the Dirac 5 represents the surface 
adhesion. Eq. (j2Jl is our working definition of the Baxter 
model. 

In the limit r — > oo, Eq. (PJ reduces to the hard sphere 
potential, so t^^ measures the strength of adhesion. Al- 
ternatively, r itself can be regarded as the effective tem- 
perature, since the true thermal energy kT has been ab- 
sorbed into the definition of the potential and does not 
appear explicitly in the right-hand side of Eq. (3). In the 
following, we shall refer to r as the temperature. 

The configuration space of an A^-particle adhesive hard 
sphere system can be notionally divided into subspaces 
within each of which the number of binary contacts is 
constant. Let 
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where Tij is the separation of particles i and j. The 
integral extends over all particle coordinates, but the 
only configurations that make a non-zero contribution are 
those in which no particles overlap and the total number 
of binary contacts is M . The latter requirement is en- 
forced by the function (5m (r^). We point out that the 
presence of delta-functions in the Boltzmann factor of a 
pair of adhesive hard spheres should not be interpreted 
as a holonomic constraint on the Lagrangian of the sys- 
tem but rather as the limit of an increasingly stiff bond. 



In fact, if the bonds between adhesive hard spheres acted 
as holonomic constraints, then the integrand in Eq. Q 
would include a Jacobian as a result of the integration 
over the (constrained) momenta (see e.g. Ref. not 
to mention the fact that the dimensionality of the phase 
space would in that case depend on the number of bonds. 
Although we will not need to manipulate an explicit form 
for 5m{y^), it may be formally expressed as a sum of M- 
fold products of 5 functions: 
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where is the separation of particles i and j . The sum is 
over all combinations of the 2M indices that give distinct 
products of 5 functions but neither generate constraints 
on a 'self-separation' (such as ru) nor duplicate a con- 
straint on any given separation within a single term in 
the sum. 

The configuration integral now takes the succinct form 
of Eq. (|SJ), which emphasizes the role of il^vM as an ef- 
fective density of states: 
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This notation will be useful when we consider histogram 
reweighting. 



III. SIMULATION 
A. Canonical Monte Carlo 

The singularity in Eq. Q means that the conventional 
Monte Carlo algorithm, in which the test particle un- 
dergoes random displacements, is not applicable. The 
probability of a trial move bringing two particles exactly 
into contact vanishes, and the energy change associated 
with breaking such a contact is formally infinite. Nev- 
ertheless, the integrals in Eq. are finite (though see 
Section [ill D|l and there is an equilibrium between states 
with different total numbers of contacts, even if the time 
scales for establishing the equilibrium diverge as the limit 
of vanishing potential range is taken^ In simulations of 
the adhesive hard sphere model it is therefore necessary 
to employ moves that explicitly make or break contacts. 

For our canonical simulations we adopt the method 
pioneered by Seaton and Glandtp*li and refined by Kra- 
nendonk and Frenkel.^" In a conventional Monte Carlo 
simulation the Boltzmann weights of configurations are 
compared, where a configuration means a specification 
of all particle coordinates. In contrast, simulations of 
the adhesive hard sphere model proceed by comparing 
the weights of different coordination states. Here we use 
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the phrase coordination state to mean a specification of 
which particles are touching each other. In general, many 
configurations are compatible with a given coordination 
state. Since the Hamiltonian depends only on the num- 
ber of contacts in the system, all configurations belonging 
to the same coordination state have the same statistical 
weight. With these concepts in mind, the key features of 
the algorithm may be outlined as follows. 

1. A particle (the "test particle" ) is chosen at random. 

2. A list is made of the possible coordination states of 
the test particle with other particles in the system, 
keeping the positions of the latter fixed. The list 
always contains one item in which the test particle 
makes no contacts at all, and — 1 items where it 
touches just one of the other particles in the sys- 
tem. It also enumerates pairs of particles that are 
sufficiently close for the test particle to touch both 
simultaneously, as well as triplets compatible with 
three simultaneous contacts to the test particle. 

3. For each coordination state in the list, a transition 
probability is calculated. 

4. A coordination state is chosen at random with 
weight in proportion to its transition probability. 

5. A uniformly distributed random configuration be- 
longing to the chosen coordination state is chosen. 

6. If the configuration generates a hard-sphere overlap 
with another particle in the system, the move is 
rejected, otherwise it is accepted. 

We now show how to calculate the transition probabil- 
ities in step 131 so that repeated application of this pro- 
cedure generates the correct Boltzmann distribution of 
configurations. The total relative unnormalized weight 
of a coordination state of the test particle is given by the 
integral of the Boltzmann factor over all configurations 
belonging to that coordination state. For the coordina- 
tion state involving no contacts with other particles, this 
weight is the free volume accessible to the test particle, 
which we label t: 



1^(0) = 




V being the volume of the simulation box. For a coor- 
dination state involving contact of the test particle with 
one other particle, a, the test particle is constrained to 
the accessible part of the surface of a: 

r ^ 

^^'^ = ^ / ^(^ta - <T)Y\eiru - a)dvt. (7) 

Similarly, for coordination states involving simultaneous 
contact with two particles (a and b) or three particles (a. 



h and c), we have 

2 r ^ 
^^'' = (li") / Kna ^ <y)5{rtb - <j)\{Q{ru - <j)dvt 
''' i^t 

(8) 

and 

W^'^ = (^^y S{rta - <y)S{rtt - <7)S{rt, - a)x 

N 

Y[Q{ru-a)drt. (9) 

Moves that create or destroy more than three contacts si- 
multaneously are not considered in the algorithm. A par- 
ticle can, however, attain a coordination number greater 
than three through suitable combinations of moves where 
three or fewer contacts are created. 

If Eqs. © to lO could be evaluated conveniently, 
the canonical distribution could be generated directly 
by choosing a coordination state at random using these 
weights. The transition probability would be 

where the possible coordination states of the test particle 
have been numbered from 1 to n, and the weight Wj of 
state j has been evaluated using the appropriate formula 
from Eqs. © to ©. Note that Uij does not depend 
on i, and that all moves would be accepted, since new 
configurations never generate overlaps and are already 
chosen with the desired probability distribution. 

The product of step functions that appears in Eqs. ® 
to ((SJ to preclude overlaps with the remaining particles in 
the system, makes evaluation of the weights very difficult. 
Following Seaton and Glandt,2ii4 we consider instead the 
modified weights, in which the overlap is initially ignored. 
The weight of the coordination state where the test par- 
ticle makes no contacts is now a free integral over the 
container volume: 

M^'(o) = f drt. (11) 
Jv 

The weight for a one-contact coordination state is an 
integral over the complete sphere of radius a surrounding 
a particle: 

^'^'^ = T?- / S{rta^a)dvt. (12) 

l^T Jv 

For two simultaneous contacts we integrate around a cir- 
cle of configurations: 

VF'(^) = (^) ' S{rta - a)S{nt - a)drt (13) 

and for three contacts, only two configurations (located 
symmetrically above and below the plane defined by a, b 
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and c) contribute: 
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(14) 

The integrals in Eqs. (|12() to H14() can be evaluated by 
suitable series of variable changes. Alternatively, one 
may consider first a square- well potential of finite width. 
Coordination states with 1, 2 and 3 contacts then restrict 
the test particle's center-of-mass to, respectively, a spher- 
ical shell, a torus with diamond-shaped cross-section, and 
a pair of parallelepipeds. Taking the limit of infinitesi- 
mal well width then produces the same results as direct 
integration of the delta functions. One obtains 



W^'(o) = V, 
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Unlike previous workerSfSii^ we consider moves to all 
triplets of particles capable of simultaneously touching 
the test particle, not just to those triplets where a, b 
and c are all mutually touching. The coordination states 
neglected by the latter approach can be reached by com- 
binations of other moves, but including the direct moves 
should make the algorithm more efficient at low tem- 
perature or high density, where the average coordination 
number is expected to be high. 

It is the weights in Eqs. I|15|) to H18|) that we use for 
the transition probability in step 01 of the Monte Carlo 
algorithm: 



W 



(19) 



To show that this choice satisfies detailed balance and 
leads to the correct Boltzmann distribution, let fj — 
Wj/Wj be the fraction of configurations belonging to co- 
ordination state j that are not blocked by overlap with 
other particles. The rejection of overlaps in step|H|of the 
algorithm means that a trial move to coordination state 
j will be accepted with probability pfj'^ — fj . The overall 
probability of moving to state j from state i is therefore 



(20) 



Note that the sum in the denominator of Eq. H2()|l de- 
pends on the positions of all the particles except the test 
particle. Since states i and j differ only by the position 
of the test particle, this sum is the same for the reverse 



move, so that Wjtt,, 



and detailed balance is 



satisfied. Recalling that we do not attempt moves that 
create more than three contacts simultaneously, detailed 



balance also requires that an attempted displacement of 
any particle that has accumulated more than three con- 
tacts is automatically rejected. 

Having chosen a coordination state in step ^ of the 
algorithm, step [S] requires a uniformly distributed ran- 
dom configuration from this state to be generated. For 
coordination states involving 0, 1, 2 or 3 contacts one 
must therefore select, respectively, a random point in the 
simulation box, a random point on a sphere surrounding 
particle a, a random point on the circle of contact with 
particles a and 6, or one of the two points lying a distance 
a from each of a, b and c. 

In practice, it is too time-consuming to enumerate all 
possible coordination states of the test particle, and this 
task scales unfavorably with the system size. It is there- 
fore desirable to restrict movement of the test particle to 
the vicinity of its original position, so that only a sub- 
set of the possible coordination states need be identified 
and considered. However, as discussed by Seaton and 
Glandt^ care must be taken since the normalizing sum 
in Eq. H2UII then depends on the original position of the 
particle. Kranendonk and Frenkel^*^ have prescribed a so- 
lution to the problem: a random point is chosen within 
a sphere of radius rtost centered on the original position 
of the test particle, and this point is used as the center 
of a second sphere (the test sphere), also of radius rtest- 
Displacements of the test particle are then considered to 
coordination states within the test sphere. The chance 
of choosing the same test sphere in the reverse move is 
identical, so that the modified scheme maintains detailed 
balance. 

The test sphere introduces another complication in 
that it may intersect some of the spherical and circu- 
lar contact surfaces available to the test particle, thereby 
excluding a fraction of the configurations belonging to 
the corresponding coordination states. One must there- 
fore calculate what fraction of each contact surface lies 
within the test sphere, and modify the weig hts in Eq. l(T^ 
by these fractions when calculating the transition prob- 
abilities in step 121 of the algorithm. In step[Sl one must 
then choose a uniformly distributed random configura- 
tion restricted to the part of the surface that lies within 
the test sphere. 



B. Grand Canonical Monte Carlo 

Since two phases can only coexist when their pressures 
are equal, computational studies of coexistence often em- 
ploy an ensemble where the size of the simulation cell can 
fluctuate in response to the pressure. The coordinates of 
the particles or the cluster centers-of-mass are typically 
scaled in proportion to the change in the cell length. In 
the adhesive hard sphere fluid, clusters can percolate even 
at moderately low density. A percolating cluster spans 
the simulation cell, so that increasing or decreasing the 
cell size in the presence of such a cluster would always 
involve breaking a contact or generating an overlap be- 
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tween at least two particles. Since overlaps are forbidden, 
and the breaking of a bond is energetically prohibited, all 
trial changes in the cell size would be rejected. Isobaric- 
isothermal and Gibbs-ensemble simulations therefore fail 
as soon as the percolation threshold is crossed. 

To study density fluctuations, we work instead in the 
grand canonical ensemble, where the volume of the sim- 
ulation is fixed, but the number of particles fluctuates in 
response to the imposed chemical potential. To avoid the 
formally infinite change in potential energy associated 
with the creation and destruction of contacts between 
particles, we only consider the insertion and removal of 
particles whose coordination number is zero. The proce- 
dure is therefore identical to that used for ordinary hard 
spheresfi^ with the slight modification that a removal 
must be rejected if the chosen particle happens to have 
a non-zero coordination number. In our grand canoni- 
cal simulations, particle insertion and removal steps are 
attempted with equal probability totaling 45%. 

In an attempt to accelerate equilibration, we also per- 
form cluster translation moves with probability 5%. A 
particle is chosen at random, and then all particles in 
the cluster to which it belongs are translated by the same 
amount. The maximum translation in each Cartesian di- 
rection is inversely proportional to the number of parti- 
cles in the cluster. 



C. Parallel Tempering 

At low temperature and high density, the formation of 
large, highly-coordinated clusters drastically slows down 
the equilibration of the simulations. To help overcome 
this problem, we employ the parallel tempering scheme 
of Geyer— in the grand canonical ensemble. In our im- 
plementation, several runs are performed simultaneously 
at the same temperature and a series of increasing chemi- 
cal potentials. At sufficiently low chemical potential, one 
can always recover an ergodic dilute gas. Parallel tem- 
pering involves periodic attempts to exchange the con- 
figurations of pairs of adjacent runs in the hierarchy of 
chemical potentials. The acceptance probability of such 
moves is given below and ensures that the correct grand 
canonical distribution is obtained for each run individu- 
ally. The advantage of swapping configurations is that 
a large cluster formed at high chemical potential breaks 
down when transferred to lower chemical potential, while 
a new cluster is built up from the configuration received 
by the run at higher chemical potential. 

Using Eq. wc can write the grand partition func- 
tion of the adhesive hard sphere model as 

CO y fj,/kT \ ^ 

N=0 ^ ^ N,M 

(21) 

where the chemical potential fi enters through the activ- 
ity z = A""^ exp(/i/fcr) and A is the thermal de Broglie 
wavelength. The probability of observing the system in 



a configuration with A'' particles and M contacts under 
conditions of effective temperature r and activity z is 

PNAiir. z) = QNMZ^T-^'/EiT, z). (22) 

The joint probability of a system having A'^i particles 
and Ml contacts at temperature r and activity zi at the 
same time as a second system with equal temperature but 
activity Z2 has N2 particles and M2 contacts is simply 
the product P = PniMi (t, zijPjss^M-i (t, 2:2)- The equilib- 
rium probability of observing a pair of systems under the 
same conditions but with the configurations exchanged 
is P' — Pn2M2{tt Zi)Ptq^Mi{T, Z2)- Since the probabihty 
of attempting a configuration exchange is independent of 
configuration, the correct joint distribution is sampled if 
exchanges are accepted with probability 

pll' - min[I, P'/P] = min[I, (zi/zz)^^"^^]. (23) 

The activities must be chosen close enough that typical 
values of p^'^ are not too small. As the extensivity of 
the exponent in Eq. (|23|l makes clear, p^'^ decreases with 
increasing system size because of the decreasing relative 
size of the fluctuations. It is therefore necessary to use a 
larger number of parallel runs with more closely spaced 
activities for larger systems. 

A natural companion of the parallel tempering tech- 
nique is multiple histogram reweighting,^^'^^ which al- 
lows the results from simulations of different conditions 
to be combined. One may then obtain rigorously inter- 
polated data for conditions different from those explicitly 
simulated without the need for further runs. Our imple- 
mentation of this method is detailed in the Appendix. 

D. Inherent Pathology 

It has been pointed out that for systems of A > 12 
adhesive hard spheres, the configuration integral Z^v 
diverges fiSiSS and that the model is therefore patholog- 
ical. The problem arises from certain maximally con- 
nected clusters in which the number of contacts between 
particles exceeds the number of vibrational degrees of 
freedom, leading to non-integrable singularities in Zjq. 

The present simulations are concerned with the fiuid 
phases of Baxter's model. Under conditions where dense 
clusters form, aggregation tends to proceed in a poly- 
tetrahedral manner, with particles attaching to existing 
triangular facets due to the contact adhesion. The clus- 
ters that lead to the divergence of are close-packed 
and characteristic of the crystal phase. In contrast, reg- 
ular tetrahedra are not space-filling, and so troublesome 
clusters are increasingly unlikely to arise spontaneously 
in the simulation as t decreases. Indeed, for such a short- 
range potential, the liquid phase is expected to be highly 
metastable with respect to the crystal^^ with a high free 
energy barrier between the two. Here we effectively ne- 
glect the diverging contributions of crystal nuclei con- 
taining 12 or more particles. 
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In the event that a pathological cluster were to arise 
in the simulation, the algorithm described in this section 
would simply fail to observe the additional contact be- 
tween particles that leads to non-integrability. Since such 
"accidental" contacts arise from the geometrical packing 
of strictly monodisperse spheres, the algorithm effectively 
introduces an infinitesimal polydispersity. For implica- 
tions regarding the applicability of the model, we refer 
readers to the detailed discussion by Stell.^^> 



IV. RESULTS 
A. Coexistence Curve 

A grand canonical simulation yields a histogram of 
the number of particles in the system — effectively the 
probability distribution of the reduced particle density, 
p = a^N/L^, where L is the length of the cubic simula- 
tion box. We work throughout with p, which is related 
to the volume fraction 77 of the spheres hy rj = p7r/6. 

At any subcritical temperature, there is an activity at 
which the probability distribution of p is bimodal with 
peaks of equal area, signifying equilibrium between a low- 
and a high-density fluid phase. The coexisting densities 
are given by the mean of each peak. In practice, the fi- 
nite size of the simulation box means that for modestly 
subcritical temperatures, the two peaks overlap, and it 
becomes impossible to attribute the intermediate densi- 
ties to one peak or the other. It is therefore more reliable 
to define equilibrium by equal heights of the peaks, and 
the coexisting densities by the peak positions. 

We note that, at significantly subcritical temperatures, 
where two separated peaks are observed, the effect of 
small changes in the activity about the coexistence point 
is to alter the relative heights of the peaks without af- 
fecting much their positions. Under these conditions, the 
positions of the peaks at coexistence are also not very 
sensitive to the system size. Closer to the critical point, 
discrepancies due to the definition of coexistence and fi- 
nite size effects both increase. However, using the peak 
heights is far less ambiguous than using their areas. 

Table Q reports the coexisting (peak) densities as a 
function of the temperature parameter r in a cubic simu- 
lation box of length L = 8a. Also given are an estimate of 
the statistical error and a measure of the size dependence 
of the results by comparison with simulations at L = 6a. 
At the lowest temperatures, where the coexisting density 
peaks are well separated, the difference between the peak 
positions and their integrated means is roughly 0.01, the 
direction of the differences being such that the means of 
the two peaks are closer than their peak positions. 

The pronounced size dependence near the top of the 
low- and high-density branches was noted in previous 
work,— and was exploited to derive the critical point 
in the thermodynamic limit, as listed in the last line of 
Table At this point, the adhesive hard sphere den- 
sity distribution can be scaled onto the well-known^'^ uni- 



TABLE I: Coexisting peak densities pio and phi of the low- 
and high-density fluid phases of the adhesive hard sphere 
model with simulation cell length L = Sct. Jran denotes the 
uncertainty in locating the peak positions. S^ys is the average 
of the decrease in pio and the increase in phi when the sim- 
ulation cell size is decreased to L = 6(j. This quantity is an 
indication of the increasing size-dependence of the results as 
the critical point is approached from below. 
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TABLE II: Critical temperature and density of the adhesive 
hard sphere model from the Percus-Yevick compressibility- 
and energy'? routes, and from simulation, 



method 




Pc 


PY compressibility 


0.0976 


0.232 


PY energy 


0.1185 


0.609 


simulation 


0.1133 ±0.0005 


0.508 ±0.01 



versal critical order-parameter distribution of the three- 
dimensional Ising model at any system size2iiS*2& (5 < 
L/a < 10 were studied). However, it is clear from Table 
^and Fig. n that the finite-system coexistence branches 
significantly overshoot the critical point. It is partly be- 
cause of the varying finite-size effects that we have re- 
frained from fitting the coexistence curve to a functional 
form, such as a scaling law with universal exponents. We 
also point out that the parameter r is not the thermo- 
dynamic temperature in the usual sense. Even if tem- 
perature is the relevant parameter in an experimental 
realization of adhesive hard spheres, its relation to r de- 
pends entirely on the substance, and could be linear(21 
inverse,^- or more complicated.— 

Table |H] and Fig. ^ emphasize the large discrepancies 
in the critical point and coexistence curve derived from 
the compressibility^ and energy^ routes of Percus-Yevick 
theory; the critical densities differ by almost a factor of 
three. The present simulations indicate that the true 
properties lie between the two theoretical results, but 
somewhat closer to the energy route, despite the fact 
that the compressibility results are more often used when 
modeling experimental data. 
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FIG. 1: Phase diagram of the adhesive hard sphere fluid. 
The two solid lines are the fluid-fluid coexistence curves from 
Percus-Yevick theory as marked. The simulation coexistence 
data from Table|I]are shown by points with error bars (mark- 
ing statistical uncertainty) , with a dotted line merely to guide 
the eye. Critical points are indicated by crosses. The dashed 
lines show the percolation threshold: long dashes are the 
Percus-Yevick result, Eq. II25II and short dashes with circles 
are the simulation results and their fit, Eq. 1241 . To the high- 
density side of the threshold the system essentially always 
contains an infinite cluster. 



B. Percolation Threshold 

In principle, a percolating cluster is a dynamical ob- 
ject. A system-spanning cluster may assemble, but its 
constituent particles are in dynamic equilibrium with the 
rest of the system. In a simulation, one can therefore 
define a percolation probability equal to the fraction of 
the encountered configurations that contain a system- 
spanning cluster. At a given temperature, this probabil- 
ity is a sigmoidal function of density, rising from zero to 
unity in a relatively narrow range. The density at which 
the probability passes through 0.5 has been found to be 
quite insensitive to the system size^S*ii and we therefore 
adopt this point as a robust definition of the percolation 
density, Ppcrc{T). For p > ppcrc the system effectively 
always contains a percolating cluster. 

We have measured Ppcrc(''") in canonical simulations 
of a system with N — 500 particles. The decreasing 
density at which percolation sets in as r is lowered means 
that the ergodicity problems associated with dense low- 
temperature regime do not arise, and parallel tempering 
is not necessary. Our data are well reproduced by a ratio 
of polynomials with six adjustable parameters: 



Ppcrc ('^) 



-10.09 + 182 At + 606. Or^ + 15.31t^ 
1 + 507.9t + 548.9t2 



(24) 



applicable in the range 0.095 < t < 3. The maximum 
residuals between the measured and fitted densities are 
0.001 for the portion above r = 0.2 and 0.002 below 
it. These errors are equal to the estimated uncertainty 
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FIG. 2: Percolation threshold of the adhesive hard sphere 
fluid. Open circles are Monte Carlo data from this work, and 
the dashed line is the fit Eq. (1241 . F illed circles are Monte 
Carlo data from Table I of Ref. lla . The solid line is the 
Percus-Yevick result. 



in the simulation data, as ascertained by independent 
simulations with different random number seeds. 

Finite-size effects are relatively small in the percolation 
curve. Reducing the system size to iV = 256 particles 
shifts the curve systematically to higher density by about 
0.001 towards the high-r end of the range studied and 
0.003 at the low-r end. 

Figure |21 shows the simulation data and the fit in 
Eq. H24|) over the full range studied. We note that, 
at high r, the threshold density enters the two-phase 
crystal-fluid region of ordinary hard spheres (0.939 < 
p < 1.038). 2^ The onset of hard sphere behavior places 
an upper limit on the value of r to which we can trace the 
percolation threshold. No meaning should be attached to 
the empirical fit in Eq. H24|) outside the specified range. 

Leeii has performed a detailed scaling analysis of per- 
colation in adhesive hard spheres in the range 0.1 < r < 
0.7, and Fig.|21shows that the latter results are consistent 
with the present work. Also plotted is the Percus-Yevick 
percolation threshold derived by Chiew and Glandt from 
the divergence of the mean cluster size with temperature 
at a given density!^ 



' pcrc 



iv) 



1977^ -2ri + l 
12(1-77)2 ■ 



•q = p7r/6. 



(25) 



This formula increasingly overestimates the percolation 
density at given temperature for r > 0.177, but under- 
estimates it at lower temperatures approaching the two- 
phase regime (Fig. . 

Importantly, Fig. shows that the critical point for 
the fluid-fluid phase transition lies well within the per- 
colated regime, in contrast to the results of the Percus- 
Yevick compressibility results that suggest the critical 
point is rather close to the threshold. In an experimen- 
tal system with very short range attractive forces, the 
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onset of infinite clusters might therefore interfere with 
equihbrium phase separation and critical bchaviorp2£ 



C. Equation of State 

In the grand canonical ensemble, the pressure can be 
derived from the average density {p[z, t)) as a function 
of the activity z and temperature t. For a change in 
chemical potential at constant volume and temperature, 
VdP = {N)dii. Integration yields an isotherm of the 
reduced pressure P* = a^P/kT as a function of the re- 
duced density: 



(p(z,T))d(lnz) 



(26) 



Below the critical point, one must integrate the mean 
densities of the two fluid phases separately, determining 
the constant of integration for the high-density branch 
by the requirement that the pressures be equal at the 
coexistence activity. 

The multiple histogram technique is invaluable in eval- 
uating Eq. H26|l . since it allows the density to be found 
effectively as a continuous function of the activity from 
a series of simulations at discrete points. Equally, one 
can interpolate between simulations at two values of r to 
obtain an intermediate isotherm. 

We have performed grand canonical parallel tempering 
simulations spanning a wide range of particle densities 
(0 < p < 1) at each of r ^ 0.1, 0.105, 0.11, 0.12, 0.13, 
0.15, 0.2, 0.25, 0.35, 0.5, 1, 1.5, 2, 3, and 5 using a sim- 
ulation box of length L — 8cr. These directly-measured 
isotherms were supplemented by 22 intermediate temper- 
atures derived by histogram reweighting from the pair of 
simulations bracketing each temperature. 

The rather unconventional definition of the adhesive 
hard sphere model in Eq. ^ has consequences for fitting 
of an empirical equation of state to the data P{p,t). 
Firstly, the high temperature limit is not an ideal gas 
but a system of hard spheres, whose equation of state is 
reproduced accurately by Speedy's formulai^ 



0.076014x2 + 0.019480^3 



1 - 0.548986X -I- 0.075647x2 



(27) 



where x = B^^p and 



2tt/3 is the reduced 
second virial coefficient of hard spheres. Additionally, 
the first four virial coefficients of adhesive hard spheres 
are known analytically, '^^ consisting of the corresponding 
temperature-independent hard-sphere part plus a contri- 
bution from the adhesion. The adhesive contributions 



Bf'\T) 



B^ 

At ' 



-5r-i 



18 



(28) 
(29) 



(-13.77358t~i -I-6.114T- 



V 4 

-1.518r-3 0.17398r-'* 

-6.33514 X 10"^"^ - 6.51271 x IQ-^rt^) 

With this information, it is already clear that many ex- 
isting empirical equations of state are not of an appro- 
priate form for this model. Even the popular modified 
Benedict-Webb-Rubin (mBWR) equatiouj?-? which em- 
ploys 33 adjustable parameters and has been successfully 
fitted to the Lennard- Jones fluid equation of statCf^ is 
unsatisfactory as it stands: the terms containing frac- 
tional powers of temperature are redundant in the light 
of Eqs. (|28|l and (|29() . and its high temperature limit is 
not flexible enough to reproduce Eq. l|77|) . 

We have attempted to adapt the mBWR equation by 
removing unnecessary terms and inserting ones more ap- 
propriate to the adhesive hard sphere model. However, 
we have been unable to obtain a satisfactory fit across 
the full gamut of our simulation data and must therefore 
be satisfied with an empirical fit to the supercritical part 
of the equation of state, which we now present. We note 
that the supercritical regime is relevant to many experi- 
mental applications of the adhesive hard sphere model. 

To guarantee the correct low-density (up to third or- 
der in p) and high-temperature limits, we write the full 
adhesive hard sphere pressure as 

P*{p, t) = P^p) + Bf'^{T)p' + Bf'\T)p^ + P(?,(p, r), 

and fit the final term to a ratio of mixed polynomials in 
p and 1/r: 



(32) 



In Eq. H32|l the sums are over the sets of terms, with 
coefficients {ay} and that have been chosen to 

provide the flexibility to reproduce the data. 

We have used a standard nonlinear least-squares 
routine^^ to optimize Eq. to reproduce our simulated 
and interpolated isotherms P* (p, t) for the L — sys- 
tem at 35 values of r > 0.117 (the apparent finite-system 
critical point lying just below t — 0.117). Starting from 
a modest number of terms, contributions were inserted 
and removed according to their ability to improve the fit, 
and the final result was a set of 35 terms, 24 in the nu- 
merator and 11 in the denominator, as listed in Tables 
mil and Hvl We emphasize that the fit is meaningless be- 
low the lowest temperature used in the parametrization, 
T = 0.117, but may be applied at arbitrarily high t since 
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TABLE III: Numerator coefficients for tfie equation of state, 
Eq. 



i 


j 




i 


j 




4 


2 


-3.787821 


6 


4 


-0.1954175 


4 


4 


-0.3184060 


6 


5 


-0.3832949 


4 


5 


0.1460575 


7 


1 


-207.9268 


4 


6 


-0.01199764 


7 


2 


301.5938 


5 


1 


-24.12376 


7 


3 


-88.02139 


5 


2 


0.6172855 


7 


4 


6.609781 


5 


4 


-1.679166 


8 


1 


114.6900 


5 


5 


0.05973859 


8 


2 


-155.9904 


5 


6 


0.01606935 


8 


3 


42.26076 


6 


1 


111.7043 


8 


4 


-4.117703 


6 


2 


-137.3363 


8 


5 


0.1512075 


6 


3 


43.98811 


8 


6 


-0.003129547 




0.0 ^ • ^ • ^ • ^ • 1 
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TABLE IV: Denominator coefficients for tfie equation of 
state, Eq. l|32|l . 



i 


j 


bij 


i 


J 




1 


2 


29.25140 


3 


2 


85.86019 


1 


3 


-3.259845 


3 


3 


-8.978804 


2 


1 


32.11475 


4 


1 


24.33100 


2 


2 


-88.09174 


4 


2 


-26.90679 


2 


3 


9.484818 


4 


3 


2.751372 


3 


1 


-57.07652 









Eq. H31I) reverts smoothly to the ordinary hard sphere 
equation of state at t = oo. 

Despite the inapphcabihty of the fit below r — 0.117, 
one can always obtain the pressure in the low density 
limit using the hard-sphere equation of state plus the 
exact adhesive contributions up to fourth order in p as 
listed in Eqs. 

P*iP:r) = P^sip) + Bf''{T)p' + Bf\r)p^ + Bf^{T)p\ 

(33) 

In fact, the low-density side of the coexistence curve is 
sufficiently rarefied that Eq. H33() is quite satisfactory for 
the full gas-like branch of subcritical isotherms, except 
for r very close to the critical point. 

Returning to the supercritical equation of state. Fig. El 
shows the pressure along the r — 0.12 isotherm to demon- 
strate the performance of the fit. Note that the simu- 
lation results can be shown as a continuous line, since 
the multiple histogram technique allows us to calculate 
the pressure on an arbitrarily fine grid. The figure also 
shows the (dimensionless) derivative {dP* /dp)r, which 
is related to the dimensionless isothermal compressibil- 
ity K* hy 1/k* = p(dP* /dp)r- The fit reproduces the 
derivative acceptably, despite the derivative not being 
used in the fitting process. In fact. Fig. O is a somewhat 
pessimistic case, since the agreement between the simu- 
lation data and fit, especially for the derivative, improves 
with increasing r. 

The actual (not fractional) difference between the sim- 
ulation and the fit is plotted for three temperatures in 



FIG. 3: Comparison of the r = 0.12 isotherm from simulation 
and as reproduced by Eqs. 13H and 1321 . Also shown is the 
dimensionless derivative [dP* / dp)T from both simulation and 
fit. In each case, the solid line is the simulation data and 
the dashed line the fit, though the two are not everywhere 
discernible. 
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FIG. 4: Discrepancy (unsealed differences) between the pres- 
sure directly from simulations and the empirical form Eq. (I32II 
for three temperatures. 

Fig. ^ By construction, the discrepancy goes smoothly 
to zero at zero density. The increasing absolute size of 
the discrepancy at high densities still represents a small 
fractional error in the light of the steeply rising pres- 
sure in that regime (see Fig. |2Jl . Approaching the critical 
point, the simulation data and fit suffer from finite-size 
effects. However, simulations with a smaller simulation 
cell of L = 6(7 produce almost indistinguishable results 
for T > 0.14, so the equation of state provided here is 
robust over a wide range of temperatures. 

We now compare the simulation isotherms with the 
predictions of Percus-Yevick theory. Expressions for 
the pressure derived from the compressibility and energy 
routes, originally due to Baxter^ and Watts et al.^ are 
given in convenient and explicit form by Barboji2& and 
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reduced density, p 



FIG. 5: Adhesive hard sphere pressure on a supercritical 
isotherm. The solid curve is the simulation result at r = 0.15, 
equivalent to t/tc — 1.3239. The dashed and dot-dashed lines 
are the Percus-Yevick compressibility and energy results, re- 
spectively, each calculated at the same t/tc using the associ- 
ated Tc- Shown by dotted lines are the ideal gas and ordinary 
hard sphere equations of state. 

Tenne.'^'' Figure [S] shows the pressure along a supercrit- 
ical isotherm at t/t^ = 1.3239, which corresponds to 
r = 0.15 using the critical temperature from simulations. 
The two Percus-Yevick results are plotted for the same 
value of the ratio t/tc, but using their respective critical 
temperatures (Table llT)l . Also plotted for comparison are 
the ideal gas equation of state and Eq. 1)27(1 for ordinary 
hard spheres. We see that the attractive component of 
the adhesive hard sphere potential can dramatically re- 
duce the pressure with respect to hard spheres even at 
high density. Indeed, the attraction at this temperature 
is strong enough to overpower the effects of excluded vol- 
ume to the extent that the adhesive hard sphere pressure 
lies below that of the ideal gas up to p w 0.9 (simulation 
data). 

Figure |H| compares subcritical isotherms at t/tc — 
0.8826, corresponding to r = 0.1 for the simulation 
data. Coexisting densities have been joined with tie 
lines. There is a gap in the compressibility equation 
curve where the expression for the pressure has no real 
solution. Scaling the temperatures according to the as- 
sociated critical points places the simulation results for 
the low-density branch pressure and the coexisting densi- 
ties between the two routes of Percus-Yevick theory. We 
note, however, that the Percus-Yevick pressure is rather 
unrealistic in the high-density branch, especially from the 
energy equation, since the pressure must rise very steeply 
at sufficiently high density. 

V. CONCLUDING REMARKS 

It is noteworthy that, despite there being at least two 
analytical equations of state for the adhesive hard sphere 




reduced density, p 



FIG. 6: Adhesive hard sphere pressure on a subcritical 
isotherm. The solid curve is the simulation result at r = 0.1, 
equivalent to t/tc = 0.8826. The dashed and dot-dashed lines 
are the Percus-Yevick compressibility and energy results, re- 
spectively, each measured at the same t/tc using the associ- 
ated Tc- Tie lines join coexisting densities. 

model, arising from the compressibility and energy routes 
of Percus-Yevick theory, it is by far more common for 
the compressibility results to be cited and used. It is 
therefore significant that, at least as far as the phase 
coexistence curve is concerned, the present results place 
the energy equation closer to the truth. 

In summary, we have presented the first firm numerical 
results for the fluid-fluid coexistence curve of the adhe- 
sive hard sphere model, as well as accurate empirical fits 
to the percolation threshold and the supercritical equa- 
tion of state. Given the popularity of Baxter's model in 
the interpretation of experimental data in colloid science, 
we hope that this information will enhance the model's 
utility. The results also constitute a reference for future 
simulations of short ranged attractive systems. 
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APPENDIX A: MULTIPLE HISTOGRAM 

Recall Eq. for the probability of observing a sys- 
tem of adhesive hard spheres with N particles and M 
contacts under specified conditions of effective temper- 
ature T and activity z. The continuous probability dis- 
tribution of the density, P{p] t, z) can be approximated 
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from the discrete probability in Eq. H22() by summing 
over M and normalizing with respect to integration over 



P 



P{p; T, ^) = y X! Pnm{t, z). 



(Al) 



M 



A full knowledge of the density of states ^nm to within 
a multiplicative constant would therefore permit calcu- 
lation of the density distribution at any combination of 
temperature and activity using Eqs. (|22|l and IjAip . 

In this appendix, we detail our implementation of 
the multiple histogram techniquei2*iS for adhesive hard 
spheres. The aim is to combine the data from several 
simulations at different temperatures and activities to 
produce an optimal estimate for the density of states 
over a wide range of N and M so that Eqs. I|22|) and 
(|A1I) can be applied. We follow closely the approach 
described by Weerasinghe and Amar^^ for the canoni- 
cal ensemble. Our specific adaptation for the adhesive 
hard sphere model and the grand canonical ensemble is 
included here for completeness. 

Consider a grand canonical simulation at temperature 
Ti and activity Zi, where the subscript i labels the simu- 
lation. One can accumulate a two-dimensional histogram 
of N and M. Let SiNM be the number of sampled con- 
figurations that had N particles and M contacts, and 



NM 



let s 

Abbreviating S(ri, Zi 
estimate 



be the total number of samples, 
to Si, the histogram furnishes the 



SiNAl/si « PNM{Ti, Zi) — 



^l-NMZi T- , 



(A2) 



where the equality is approximate because of the statisti- 
cal uncertainty inherent in the simulation. If ribins com- 
binations of N and M are sampled, the corresponding 
ribins values of ^MAi/'^i can be obtained from Eq. (|A2p 
by straightforward division. The points lying in the tails 
of the distribution suffer the largest uncertainty due to 
poor sampling. 

If TT-sim simulations are performed at different condi- 
tions (ri,Zi), 1 < i < Usim, producing risim overlapping 
two-dimensional histograms, we obtain rtbins x ^sim in- 
stances of Eq. l|A2p . The task of the multiple histogram 
technique is therefore to find the ngim values of Sj and 
the ribins values of ^Inm that minimize the overall dis- 
crepancy between the probabilities they define through 
Eq. (|22|) and those measured in the various simulations. 
By treating these nsim + ?^bins unknowns as individu- 
ally adjustable parameters, we will obtain the density 
of states and partition function without making any as- 
sumptions about their functional forms. 

We start by taking the logarithm of Eq. IIA2p : 

In S'ijvM — In Si w Ih^Inm + N In Zi — M In — In . (A3) 



The left- and right-hand sides of Eq. (|A3p are the mea- 
sured and fitted values of In Pnm , respectively. The dis- 
crepancy is the residual, 



NM 



Inp, 



NM 



NM 



X, 



NM 



InS,), (A4) 



where we have introduced the abbreviations pinm = 
SiNM/si and XijvM = NliiZi- Mlnxj. 

If we assume that sampled configurations are uncor- 
related, the variance of SiMM measured in a simulation 
of finite length is expected to bei^ cr'^{SiMM) = {SiNM), 
where the angle brackets indicate the true ensemble av- 
erage over simulations of the same length. If we ap- 
proximate this average by the value actually measured, 
then we can estimate the standard deviation of SiNM 



to be cf{SiNM) = ^iNM- Propagating errors, we find 
a (In piNM) — Sifqf^i- Summing the squares of the resid- 
uals divided by their standard deviations, we arrive at 
the maximum likelihood estimator 



EE 

i N,M 



SiNMRiNM ■ 



(A5) 



The first sum in Eq. (|A5|) runs over the Usim simulations, 
and the second runs over all bins in the two-dimensional 
histogram. 

We now treat the Usim values of and the nbins values 
of ^NM as adjustable parameters to be fitted by mini- 
mizing x^- Differentiating Eq. (|A5p with respect to the 
fitting variables and equating to zero yields 

dxVdlnnNM = ~2Y,S^NMR^NM= 0, (A6) 

i 

dxVdlnE, = 2 S,NMR^NM = 0. (A7) 

N,M 



Substituting Eq. (|A4|) into Eq. IjAGp and solving for 
Inf^TVA/, we obtain Tibins equations of the form 

In VtNM = ^ qiNM {'^'D-PiNM - X^NM + In Si) , (A8) 



where qiNM = S^nm/ J2j SjNM- Similarly, substituting 
Eq. IIA4II into Eq. ljA7|l and solving for In Si, we obtain 
JT-sim equations of the form 

In Si = ^ KJVA/(ln ^Inm - In PijvM + XiNM)- (A9) 

n,m 

We now have expressions for the partition functions 
exclusively in terms of the density of states and vice 
versa. Since there are fewer of the former than the lat- 
ter, we solve for the partition functions first. Substituting 
Eq. HA8|I into Eq. (|X9|) and yields 



lnS.= ^ 



PiNM 



N,M 



— \npiNM + XiNM 



QjNMi^T^PjNM — XjNM + InSj) 



(AlO) 
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Gathering terms in on the left-hand side gives 
In Si - ^ PiNM ^ qjNM In Sj 

N,M j 

= ^ PiNM ''^{QjNM — Sij){lrLpji^M — Xji^m), 
N.M j 

(All) 

where Sij is the Kronecker deha. 

Eq. l|Alip is a set of risini hnear equations in risim un- 
knowns, and can be written in matrix notation as 

AY = B, (A12) 

where Yj — In S,; . The elements of the matrix A are 



^ij — Sij — ^ PiNAlQj 



NM, 



(A13) 



N,M 



and the elements of the vector B are 



Bi — ^ PiNM ''^{ijNM ~ %)(lnpj7VM — AjTVA/)- 

N,M j 

(A14) 



One practical point in evaluating A and B is that many 
of the SiNM are zero, making Inp^jvAf and Ing^jvAf in- 
calculable. However, such logarithms always appear in 
products involving SiNM itself. Since Imix-^o x\nx — 0, 
bins with zero contents simply make no contribution to 
Eqs. (EH), PT3|l and (|XT4jl . 



Solution of Eq. (|A12|I leads to numerical values for the 
partition functions, which can then be substituted into 
Eq. (jA8|) to obtain the density of states. In principle, 
however, one cannot obtain the absolute values of the 
partition functions from a set of equations like (|22|l . only 
their ratios. For this reason, the matrix A is singular, 
and Eq. (|A12p should be tackled using singular-value 
decomposition;^^ yielding S,; and i^NM up to a multi- 
plicative constant. 
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